Wang et al. BMC Genomics 201 3, 14:1 1 2 
http://www.biomedcentral.eom/1 471 -21 64/1 4/1 1 2 



^BMC 



Genomics 



RESEARCH ARTICLE Open Access 



Cross-species transcriptional network analysis 
reveals conservation and variation in response to 
metal stress in cyanobacteria 

Jiangxin Wang 1,2 , Gang Wu 3 , Lei Chen 1,2 and Weiwen Zhang 1,2 * 



Abstract 

Background: As one of the most dominant bacterial groups on Earth, cyanobacteria play a pivotal role in the 
global carbon cycling and the Earth atmosphere composition. Understanding their molecular responses to 
environmental perturbations has important scientific and environmental values. Since important biological 
processes or networks are often evolutionarily conserved, the cross-species transcriptional network analysis offers a 
useful strategy to decipher conserved and species-specific transcriptional mechanisms that cells utilize to deal with 
various biotic and abiotic disturbances, and it will eventually lead to a better understanding of associated 
adaptation and regulatory networks. 

Results: In this study, the Weighted Gene Co-expression Network Analysis (WGCNA) approach was used to 
establish transcriptional networks for four important cyanobacteria species under metal stress, including iron 
depletion and high copper conditions. Cross-species network comparison led to discovery of several core response 
modules and genes possibly essential to metal stress, as well as species-specific hub genes for metal stresses in 
different cyanobacteria species, shedding light on survival strategies of cyanobacteria responding to different 
environmental perturbations. 

Conclusions: The WGCNA analysis demonstrated that the application of cross-species transcriptional network 
analysis will lead to novel insights to molecular response to environmental changes which will otherwise not be 
achieved by analyzing data from a single species. 
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Background 

Many biological systems operate in a similar manner 
across a large number of species or conditions [1]. 
Cross-species analysis of genomic sequences has made 
fundamental contributions to modern biology in defin- 
ing the putative function of new genes [2]. Recent expo- 
nential growth in microarray expression datasets allows 
researchers to combine expression experiments from 
multiple species to identify genes that are not only 
conserved in sequence but also operated in a similar way 
[3]. In contrast to the static sequence measurements, 
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microarrays measure the dynamic, condition-specific 
responses of complex biological systems [4-7], and their 
comparative analysis has led to improvements in anno- 
tating gene function and inferring the evolution of inter- 
action and regulatory networks [3]. In addition, cross- 
species microarray analysis also allowed identification of 
core transcriptional networks, and their conservation 
and variation in closely related species [8,9]. 

Correlation networks are increasingly being used in 
various bioinformatics applications [10-12]. Among 
them, Weighted Gene Co-expression Network Analysis 
(WGCNA) is a computational method to describe the 
correlation patterns among genes across transcriptomic 
datasets [11,13]. WGCNA has been used in identifying 
functional clusters (modules) of highly correlated genes, 
summarizing such clusters using the module eigengene 
or an intramodular hub gene, relating modules to one 
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another and to external sample traits (using eigengene 
network methodology), calculating module membership 
measures in many systems, and correlation network ana- 
lysis can also be used to determine biological correlated 
candidate biomarkers or disease therapeutic targets [11]. 
The WGCNA approach has been applied in many 
eukaryotic systems, successfully linking molecular 
targets to oncogenic signals [14], complex traits [15], 
analyzing network divergence between human and chim- 
panzee neural patterns [16], and even comparing cross- 
species gene expression in animals recently [17]. A re- 
cent comparative study of different network analysis 
methods indicated that WGCNA could be used not only 
for constructing gene networks, but also for detecting 
modules/sub-networks, identifying hub genes and se- 
lecting candidate genes as biomarkers, using Escherichia 
coli as an empirical sample [18]. 

As important primary producers and significant 
contributors of fixed carbon budget in many terrestrial 
and marine environments, cyanobacteria have been 
present in many different environments from coast to 
open sea for -2.5 billion years [19]. Currently the marine 
cyanobacteria, Synechococcus and Prochlorococcus spe- 
cies are together responsible for at least 20% of global 
carbon fixation [20]. They play such significant roles in 
the carbon cycle that it is essential to understand to 
what environmental stresses they are susceptible and 
how they respond. Metal is one important environmen- 
tal factor for cyanobacteria. For instance, iron is required 
for photosynthesis, and its limitation and restriction of 
primary productivity have been reported and has been 
considered as one important factor in the ecology of 
cyanobacteria [21,22]. Early studies have found the limi- 
tation of Prochlorococcus cell division rates by iron in 
the equatorial Pacific sea [23], and variations in the 
abundance of Prochlorococcus iron-related genes be- 
tween oceans [24]. A recent study on transcriptomic re- 
sponse of high- and low-light-adapted Prochlorococcus 
strains to changing iron availability suggested high diver- 
gence of gene expression within strains [25]. The physi- 
ology and biochemistry associated with the iron-limited 
continuous culture of the halotolerant cyanobacterium 
Synechococcus PCC 7002 were also examined [26], and 
the regulatory network for acclimation of the obligate 
photoautotrophic fresh water cyanobacterium S. 
elongatus PCC 7942 to iron limitation was studied by 
transcript profiling recently [27]. Interestingly, oceanic 
Synechococcus strains were found much more sensitive 
to iron limitation than coastal strains [28]. Copper is an- 
other metal commonly present in natural environments. 
It has been reported that different cyanobacteria species 
have varying levels of copper tolerance, with 
Prochlorococcus species considered to be more sensitive 
to copper than Synechococcus species [29]. Within 



Synechococcus species, copper tolerance has also varied 
significantly [30], with coastal species of Synechococcus 
species exhibiting increased tolerance to copper shock 
and a distinctive transcriptional responses relative to 
those of open-ocean species [31]. 

To seek better understanding of molecular response of 
cyanobacteria to metal stress, as a proof of concept we 
applied a cross-species transcriptional network analysis 
to four important cyanobacteria species under two dif- 
ferent conditions of metal stress, Prochlorococcus MED4 
(PMM) and Prochlorococcus MIT9313 (PMT) under 
stress of iron depletion [25], Synechococcus 9311 (SYG) 
and Synechococcus WH 8102 (SYW) under stress of high 
copper concentration [31]. Transcriptional networks for 
each individual species were constructed using a 
WGCNA method [11,13], and then cross-strain analysis 
was performed to reveal the conservation and variation 
in terms of the response to iron and copper stresses 
among all species. In this study, cross-species analyses of 
two Synechococcus species under high copper stress con- 
dition, and two Prochlorococcus species under iron de- 
pletion stress condition were first performed separately, 
the results revealed the signature responses to iron in 
Synechococcus and to copper in Prochlorococcus, respect- 
ively. Then a cross-species comparison was conducted 
for all four cyanobacteria species under two different 
types of metal stresses. Notably, the results showed that 
9 genes were commonly regulated in all four species 
used in this study. Although still needs more experimen- 
tal evidences, the set of genes may represent an import- 
ant core signature response to metals in these four 
cyanobacterial species. Interestingly, the species-specific 
hub genes detected in four species showed no overlap 
between each other, indicating possible species-specific 
strategy for metal acclimation in each species. Thus, the 
analysis demonstrated that the application of cross - 
species transcriptional network analysis could lead to 
novel insights into molecular response to environmental 
changes which will otherwise not be achieved by analyz- 
ing data from a single species. 

Results and discussion 

Construction of cyanobacteria metal response networks 

Comparison of transcriptional networks between differ- 
ent cyanobacteria species could provide valuable insights 
into the core and important responses to environmental 
stress at the molecular levels. To this end, we first 
sought to compile a complete transcriptional profiling 
datasets of different cyanobacteria species under iron de- 
pletion and copper toxicity. The datasets used in our 
analysis included transcriptional measurements of a total 
of 50 samples from 4 cyanobacteria species (i.e. PMM, 
PMT, SYG, SYW), representing iron depletion (PMT 
and PMM) and high copper toxicity (SYG and SYW) 
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treatments, respectively (Table 1) [25,31]. To ensure that 
our dataset were reasonably selected, we also compared 
the treatment sample datasets with their controls 
through a clustering analysis. The results showed that all 
stress-treated datasets tended to be grouped together 
and can be visibly separated from the controls, 
suggesting the stresses of iron depletion and copper tox- 
icity have caused significant changes at the transcrip- 
tional level (Additional file 1: Figure SI). Next, we 
applied the datasets of each species for the transcrip- 
tional network construction using a WGCNA method. 
Analysis of the transcriptional networks showed that a 
total of 17, 21, 32 and 28 distinct transcriptional 
modules can be detected within the transcriptional 
networks of PMT, PMM, SYG and SYW, respectively 
(Additional file 2: Figure S2). More transcriptional 
modules detected in SYG and SYW than PMT and 
PMM were probably due to the relatively large genome 
sizes of Synechococcus [32-34]. The association analysis 
between the phenotypes {i.e. iron or copper stress) and 
the detected modules showed that the distinguished 
transcriptional modules highly associated with the 
phenotypes can be identified in each of the four species 
according to the correlation coefficients (r value) and 
their confidence (low ^-values) (Additional file 3: Figure 
S3). Among all transcriptional modules detected in each 
species, we found that 5 of 17 and 6 of 21 detected 
modules significantly correlated with iron depletion con- 
dition in PMT and PMM, respectively (Table 2); and 12 
of 32 and 11 of 28 detected modules significantly 
correlated with the high copper toxicity in SYG and 
SYW, respectively (Table 2). The KEGG pathway ana- 
lysis of these phenotype-correlated modules in different 
species showed that a total of 62-79 KEGG pathways 
were involved in these modules (Table 2). Interestingly, 
69 and 49 KEGG shared pathways were detected in 
PMT/PMM to iron depletion and SYG/SYW to high 
copper toxicity, respectively (Table 2). 

Several cyanobacteria species produce secondary 
metabolites such as xanthophyll carotenoids and their 
synthetic genes are often up-regulated in response to 
high light intensity, ultraviolet radiation, and desiccation 
[35]. The ability to exist in two redox states makes iron 
an essential cofactor for proteins involved in numerous 
major cellular processes such as respiration, amino acid 
metabolism and DNA metabolism. For instance, iron 



starvation down-regulated most amino acid biosynthesis 
related genes since many enzymes are iron dependent 
[36,37]. It was also reported that iron starvation also 
induced DNA recombination and DNA repair process in 
a Gram-negative diplococci bacteria [38]. PMT was ori- 
ginally isolated from the gulfstream 135 meter in depth 
(low light adapted) and could grow with iron of an order 
of magnitude lower [25]. Analysis of the iron-responsive 
modules suggested that the significant down-regulation 
of "amino acid metabolism" and "microbial metabolism 
in diverse environments" probably reduced the iron re- 
quirement, and the up-regulation of the "metabolism of 
secondary metabolites" protected the cells from oxida- 
tive stress caused by iron depletion directly or indirectly. 
At the same time, the up-regulation of "homologous re- 
combination" might result in an increase of DNA mis- 
match repair and recombination [39]. Detailed analysis 
of genes with the same change trends (up- or down- 
regulated) in both PMM and PMT under iron depletion 
condition suggested the mechanism that cyanobacteria 
employ to iron deficiency could be: firstly transporters 
were activated to obtain more iron from the environ- 
ment [25], and photosynthesis and chlorophyll meta- 
bolism was inhibited due to the limited supply of iron 
[40-42], and then oxidative stress was enhanced and 
may cause RNA degradation, protein translation inhibi- 
tion and other cellular metabolic changes. All these fur- 
ther resulted in inhibition of DNA replication and cell 
cycle block, and eventually growth cease. 

The copper-responsive SYW yellow module showed a 
significantly high correlation between "ribosome biosyn- 
thesis and assembly" (syw03010, total 20 related genes 
enriched) and the copper toxicity treatments (p = 6.78E- 
15), suggesting an involvement of protein synthesis pro- 
gram in copper acclimation (Additional file 4: Table SI). 
Another copper-responsive SYW greenyellow module 
showed a significant correlation with the "aminoacyl- 
tRNA biosynthesis" and "citrate cycle" (TCA cycle) 
pathways, which was characterized as a down-regulation 
response. According to KEGG pathway analysis, genes in 
this module tended to encode proteins involved in 
amino acid metabolism (p = 0.0012-0.0372) (Additional 
file 5: Table S2). The genes involved in "polycyclic aro- 
matic hydrocarbon degradation" were identified in the 
copper-responsive SYG darkgreen module. The module 
also showed a significant correlation with amino acid 



Table 1 Cyanobacteria species and related datasets used in this study 


Species 


Code 


Condition (GEO ID) 


Original sample size 


Included in this study after filtering 


Synechococcus 931 1 


SYG 


Copper toxicity (GSE13910) 


14 


14 


Synechococcus WH 8102 


SYW 


Copper toxicity (GSE13910) 


10 


10 


Prochlorococcus MED4 


PMM 


Iron depletion (GSE26533) 


16 


10 


Prochlorococcus MIT 9313 


PMT 


Iron depletion (GSE26533) 


23 


16 
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Table 2 Modules and KEGG pathways significantly 
responsive to the metal treatments 



Metal 
stress 


Species 


Modules related to 
phenotype 


KEGG 
Pathways 


Common 
pathways 


Iron 


PMT 


5 


79 


65 




PMM 


6 


68 




Copper 


SYW 


12 


63 


49 




SYG 


11 


62 





metabolism (histidine, tyrosine, cysteine and methionine, 
arginine and proline; p = 5.59E-05-0.005). KEGG path- 
way analysis indicated that the SYG darkgreen module 
was enriched with genes encoding photosynthesis an- 
tenna proteins, and genes related to "microbial metabol- 
ism in diverse environments" and "oxidative 
phosphorylation" (Additional file 6: Table S3). The 
copper-responsive SYG brown module showed signifi- 
cant correlation with the "curated two-component", 
overrepresented with genes participating in signal trans- 
duction activity, including "response regulator", "two 
component system", "histidine kinase" and "ABC 
transporters". (Additional file 7: Table S4), suggesting 
the strong involvement of signal transduction systems in 
copper response in SYG species. 

For high copper toxicity, a recent study showed that 
costal strains of marine Synechococcus species (SYG in 
this study) exhibited increased tolerance to copper shock 
and a distinctive transcriptional response relative to 
those of open-ocean strains (SYW) [31]. SYG was 
predicated to be more metal-tolerant and may be better 
in adapting to fast changing environments, including bi- 
otic and abiotic environmental stresses besides copper. 
Early studies showed that toxicity of excess copper may 
be arose through several mechanisms: production of re- 
active oxygen species (ROS), metal competitive binding 
and photosynthesis inhibition [31,43,44]. In this study, 
we found that significant copper-responsive modules 
existed in both species of SYG and SYW. Modules of 
SYW greenyellow and SYG darkgreen both contained 
genes related to "amino acid metabolism" and "microbial 
metabolism in diverse environments", representing a 
general stress response. Meanwhile, more specific 
responses were also found in SYW yellow module which 
contained a down- regulated cluster of ribosomal genes, 
and in SYG brown module which had up-regulated sig- 
nal transduction and response, such as "curated two- 
component system", RR and HK system. The overlap of 
orthologous genes between SYW and SYG was modest, 
with only 46 orthologous genes being significantly 
regulated by high copper. Based on these results, we 
proposed that a possible strategy for SYW cells to sur- 
vive under high copper condition was to slow down de 
novo protein synthesis and employ oxidative protection 



to eliminate the ROS damage caused by copper treat- 
ment, while SYW cells may use relatively active signal 
transduction and response systems to adjust its metabol- 
ism to adapt the changing environments. 

Signature genes and pathways of different cyanobacteria 
species under metal stresses 

Using WGCNA, we obtained modules that contained an 
exact number of assigned genes. The same genes could 
be assigned to multiple modules, although the associ- 
ation strengths in different modules could vary. Our 
results showed a high degree of cross-species module 
similarity between PMM and PMT under iron depletion 
condition, and between SYW and SYG under high cop- 
per condition, suggesting possible core responses to each 
treatment. In addition, the overlaps between iron deple- 
tion and copper toxicity response modules were also 
observed, suggesting the possible core and general res- 
ponse to metals and possible cross-talking response net- 
works between these four species to different metal 
stresses. We identified the shared genes between response 
modules in different cyanobacteria species based on the 
significance of module membership values (Additional 
file 8: Figure S4, Additional file 9: Table S5 and Additional 
file 10: Table S6). The results showed that a total of 440 
and 430 genes were found correlated with iron depletion 
in PMM and PMT, respectively, with only 34 of them 
shared between PMM and PMT with confident low 
^-values (Additional file 8: Figure S4A). Similarly, the 
numbers of the copper correlated genes were 736 and 
580 in SYG and SYW, respectively (Additional file 8: 
Figure S4B), with only 70 of them shared between SYG 
and SYW. 

In a recent study, only 4 of 1159 orthologous genes of 
PMT and PMM were found differentially regulated in re- 
sponse to iron in both species (using fold change > 2.0), 
although expression level of over a hundred genes 
changed in two species, suggesting a great diversity in 
terms of iron adaptation among different Prochlorococcus 
species [25]. Interestingly, only 1 pair (petF, ferredoxin, 
PMM0898/PMT1429) of 4 differentially expressed genes 
(petF, isiB PMM1171/PMT0801, hliOS/08 PMM1404/ 
PMT1154, and idiA PMM1164/PMT0287) in iron stress 
previously reported [25] was showed in our shared gene 
list, due to the low correlation coefficient of the other 
three genes in the modules with r ranging from 0.001 
to 0.42. Conventionally for microarray analysis, a cutoff 
of > 2-fold change was commonly used; it is quite possible 
that some of the genes with important functions will 
have a smaller change (fold change < 2.0) and thus 
will be excluded. In this study, we focused on the co- 
expression property across two different species, thus 
no fold change threshold was applied, by doing this 
we could avoid information loss due to artificial fold- 
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change cutoff (i.e. genes changed with significantly low 
p values, but with less than 2.0 fold changes). Instead, 
we focused on only those genes in highly correlated 
modules to explore their possible functions in metal 
acclimation. 

Interestingly, our cross-species analysis revealed 9 com- 
monly shared genes among these four cyanobacteria 



species under two different metal stress conditions (Table 3). 
They were ABC transporter ATP binding component, 
carboxypeptidase Taq metallopeptidase, cell division protein 
(ftsW), ferredixin, isocitrate dehydrogenase (NAD + ) {idh3)> 
LysM domain, phosphoribosylaminoimidazolecarboxamide 
formyltransferase/IMP cyclohydrolase (purH), phycocyanob 
ilinferredoxin oxidoreductase (pcyA) and 6-pyruvoyl 



Table 3 Genes responsive in all four cyanobacteria species 



KOJD 


Name 


Taxon 


Gene 


Module 


Module 
cor 


Fold change 
(log2) 


p-value 


K01737 


6-pyruvoyl tetrahydrobiopterin synthase 


PMM 


PMM0106 


turquoise 


-0.755 


-0.5 


0.024 






PMT 


PMT0164 


turquoise 


-0.61 1 


0.3 


0.0009 






SYW 


SYNW1504 


black 


0.791 


1.5 


0.013 






SYG 


sync_1898 


brown 


0.977 


0.2 


0.002 


K00602 


phosphoribosylaminoimidazolecarboxamide formyltransferase / 


PMM 


PMM0266 


magenta 


0.714 


0.3 


0.028 




IMP cyclohydrolase 


PMT 


PMT1857 


purple 


0.613 


0.3 


0.025 






SYG 


sync_0289 


turquoise 


-0.806 


-0.4 


0.0001 






SYW 


SYNW0249 


greenyellow 


-0.957 


-0.6 


0.029 


K00030 


isocitrate dehydrogenase (NAD+) 


PMM 


PMM1596 


magenta 


0.714 


0.4 


0.008 






PMT 


PMT1935 


turquoise 


-0.61 1 


-0.2 


0.006 






SYG 


sync_0214 


pink 


0.757 


-0.5 


0.002 






SYW 


SYNW0166 


yellow 


-0.964 


-0.5 


0.001 


K03588 


cell division protein FtsW 


PMM 


PMM1458 


green 


0.661 


-0.6 


0.023 






PMT 


PMT1475 


turquoise 


-0.61 1 


-0.2 


0.027 






SYG 


sync_2326 


turquoise 


-0.806 


-0.3 


0.012 






SYW 


SYNW0475 


yellow 


-0.964 


-0.2 


0.011 


K05371 


phycocyanobilinferredoxin oxidoreductase 


PMM 


PMM0747 


cyan 


0.743 


-0.3 


0.0003 






PMT 


PMT0590 


pink 


0.695 


0.3 


0.026 






SYG 


sync_1656 


turquoise 


-0.806 


-0.4 


0.001 1 






SYW 


SYNW1084 


brown 


0.869 


0.2 


0.048618 


K01299 


carboxypeptidase Taq (M32) metallopeptidase 


PMM 


PMM0493 


red 


-0.813 


-0.8 


0.015 






PMT 


PMT1279 


greenyellow 


-0.955 


-0.3 


0.0005 






SYG 


sync_2040 


purple 


-0.808 


-0.4 


0.0008 






SYW 


SYNW1098 


greenyellow 


-0.957 


0.3 


6.78E-05 




LysM domain 


PMM 


PMM0330 


turquoise 


-0.755 


0.7 


0.003 






PMT 


PMT0190 


turquoise 


-0.61 1 


-0.2 


0.004 






SYG 


sync_0565 


brown 


0.977 


-0.7 


2.19E-05 






SYW 


SYNW1957 


yellow 


-0.964 


-0.3 


0.0007 




ABC transporter, ATP binding component 


PMM 


PMM0290 


red 


-0.813 


0.4 


0.002 






PMT 


PMT1635 


turquoise 


-0.61 1 


-0.3 


0.005 






SYG 


sync_2588 


darkturquoise 


-0.746 


-0.2 


0.0207 






SYW 


SYNW0320 


black 


0.791 


0.7 


0.015 


K02639 


ferredoxin 


PMM 


PMM0898 


red 


-0.813 


-0.1 


0.035 






PMT 


PMT1429 


purple 


0.613 


0.6 


0.011 






SYG 


sync_1953 


turquoise 


-0.806 


-0.6 


0.006 






SYW 


SYNW1277 


midnightblue 


-0.722 


0.4 


0.012 
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tetrahydrobiopterin synthase (pts) (Table 3). However, not 
all of these genes shared the same expression patterns cross 
species, except for one gene, ftsW, which was down- 
regulated in all four species under both iron starvation and 
copper toxicity. Encoding a cell-division intergral mem- 
brane protein, ftsW has been found essential in E. coli [45] 
and was required for Z-ring stabilization during sporula- 
tion septation in Streptomyces coelicolor [46]. Both FtsQ 
(another cell division protein) and FtsW are indispensable 
to Synechocystis and their depletion led to slow growth 
and giant cells [47]. Recent study indicated that FtsW was 
also directly involved in DNA damage checkpoint 
coordinately interacting with an SOS regulon in 
Caulobacter crescentus [48], and the down- regulation of 
this gene might inhibit cell division [49]. Both iron deple- 
tion and high copper can cause oxidative stress and DNA 
damage [25,31], it could thus be speculated that ftzW gene 
may play an important role in cellular defense to environ- 
mental disturbances. 

Several other stress-regulated genes detected in all four 
selected cyanobacteria species also showed different di- 
rection of regulation {i.e. up- or down-regulation) in all 
species. For instance, ABC transporter ATP binding com- 
ponent [50] was up-regulated in PMM and SYW but 
down-regulated in PMT and SYG; carboxypeptidase Taq 
metallopeptidase was up-regulated in SYW but down- 
regulated in the other 3 species. Isocitrate dehydrogenase 
(NAD + ) (idh3), involved in several metabolism pathways 
like citrate cycle and secondary metabolite biosynthesis, 
was up-regulated in PMM but down- regulated in PMT, 
SYG and SYW. Similar to idh3, LysM domain gene en- 
coding a LysM-containing protein involved in bacterial 
cell wall degradation [51], was also up-regulated in PMM 
but down-regulated in PMT, SYG and SYW. Interestingly, 
gene encoding PurH catalyzing the last two steps of de 
novo purine biosynthesis in bacteria and eukaryotes [52], 
was up-regulated in iron depletion but down-regulated in 
copper toxicity treatments. 

Phycocyanobilimferredoxin oxidoreductase, PcyA, ca- 
talyzes biosynthesis of the phycobili-protein and phyto- 
chrome chromophore precursor phycocyanobilin [53]. 
Phycocyanobilin can serve as a light-harvesting pigment 
in the photosynthetic light-harvesting structures of cyano- 
bacteria called phycobilisomes [54]. Its encoding gene was 
found up-regulated in PMT and SYW but down-regulated 
in PMM and SYG. 

Ferredoxin coding gene was up-regulated in PMT and 
SYW but down-regulated in PMM and SYG. Level of 
cellular ferredoxin and flavodoxin was previously pro- 
posed as an indicator of iron stress in cyanobacteria 
[55]; however, a recent analysis about iron availability and 
production of ferredoxin in Antarctic sea ice diatoms 
showed no strong correlation between their levels and 
stress conditions [56]. Our results here supported the 



argument that ferredoxin level alone may not be a strong 
evidence of iron limited growth. 

Another gene responsive in both iron and copper 
perturbations was pts encoding an enzyme involved 
in the early biosynthetic pathway of pteridines in 
cyanobacteria [57]. The gene was found up-regulated 
in species of PMT, SYW and SYG. The role of cya- 
nopterin in UV/blue light single transduction of 
cyanobacterium Synechocystis sp. PCC 6803 photo- 
taxis was verified recently [58]. Cyanobacteria are 
well known for producing high amounts of pteridine 
glycosides, and pterin compounds have been sugges- 
ted as possible photoreceptor pigments in some stress 
responses that are induced by UV and blue light, 
probably due to their chemical and photophysical pro- 
perties [59]. 

With the aid of WCGNA network analysis, it is also 
possible to obtain novel insights into molecular respon- 
ses to environmental changes which will otherwise not 
be achieved by analyzing data from a single species. For 
example, a total of 27 shared metal-responsive KEGG 
pathways between two metal stresses were determined 
by the WGCNA analysis (Additional file 11: Table S7), 
which included 20 pathways of "amino acid, nitrogen, 
sulfur metabolism", 2 pathways of "photosynthesis and 
porphyrin and chlorophyll metabolism", 1 pathway each 
for "ABC transporters", "oxidative phosphorylation", 
"DNA mismatch repair", and "biosynthesis of secondary 
metabolites", respectively. Based on the fact that they 
responded to both iron depletion and high copper tox- 
icity in multiple cyanobacteria species, we proposed that 
these pathways may represent the core metal-responsive 
pathways. Moreover, among all 40 signature genes of the 
PMT-greenyellow module, many were involved in amino 
acid metabolism, secondary metabolites, ABC transpor- 
ters, and homologous recombination. PMM, originally 
isolated from surface water of the Mediterranean ocean, 
is high-light adapted with higher copper tolerance [25]. 
Probably due to these physiological differences, some 
highly correlated modules detected in PMT were not 
seen in the marine coastal species PMM, although a 
large number of active KEGG pathways detected were 
commonly shared by two species under the iron deple- 
tion condition (Table 2). It has been speculated that the 
better adaptation of PMT to lower iron concentration 
may be due to more efficient iron transport systems, bet- 
ter protection of iron stress and lower requirement for 
iron during growth [25]. 

Divergent hub genes for different cyanobacteria species 
under metal stress 

Hub genes which have high connectivity within the 
networks are considered to play important roles in dif- 
ferent phenotypes, such as yeast viability [60] and 
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drought stress [61]. Recent work on cyanobacteria 
showed that there existed big divergence of gene expres- 
sion patterns between cyanobacteria species from differ- 
ent ecological environments [25,31]. 

The genes might be considered as hub genes to a spe- 
cific stress with the most connected genes in the respon- 
sive modules developed from the WGCNA analysis 
[60,61]. One of the major goals of this cross-species ana- 
lysis was to discover the hub genes in each cyanobacteria 
species under specific environmental stimuli. The 
WGCNA analysis showed that some genes were 
connected with more than 8 other genes. Based on the 
annotated genes in the WGCNA modules, we obtained 
modest lists of orthologous hub genes for PMT respon- 
sive to iron depletion (number of the genes: 16) and 
SYW to copper toxicity (16), and shorter lists for PMM 
(10) and SYG (6). Surprisingly, no overlap of the hub 
genes was found between species, suggesting they may 
be unique to each species. For example, in PMT, the 
hub genes mainly encode metabolic enzymes, like GMP 
synthase (glutamine-hydrolysing), ketol-acid reductoi- 
somerase; while in PMM, the genes were involved in 
photosynthesis system (psal and rbcL) and two-com- 
ponent system (phoB). Nearly half of the hub genes 
in SYW encode ribosomal proteins (7) and oxidative 



phosphoraylation proteins (atpfla and atpfld), com- 
pared to SYG hub genes with functions of metabolic 
enzymes and photosysm II PsbK protein and phyco- 
erythrin-associated linker proteins (Additional file 12: 
Table S8). In terms of the functionally unknown genes in 
the modules (Additional file 13: Table S9), the analysis 
showed that they were also included in some gene 
clusters. For instance, PMT101-102-107 (black module), 
PMM1395-1396 (green module), sync_2434-2436 (pur- 
ple module) were grouped together, while SYWN0193- 
0194 and SYWN0319-0320 were not in the same 
module. The function of these stress-responsive gene 
clusters with unknown functions may worth further 
investigation [60]. 

Hub genes obtained through the WGCNA analysis 
also showed distinct genes involved in iron stress in dif- 
ferent Prochlorococcus species, which provided import- 
ant clues for the differential fitness of these two species 
at low iron concentrations: i) the existence of a unique 
gene, lipopolysaccharide transport system ATP-binding 
protein (ABC-2.LPSE.A) and the activation of a urea 
transport system permease protein (urtB) may cause a 
more efficient iron transport features in PMT; ii) down- 
regulated ribosome genes (rp-S15, rp-L3) in PMT could 
save energy for cell survival under low iron 



• Iron 




Figure 1 Schematic representation of responsive genes regulated by metals in cyanobacteria. See text for details, based on WGCNA 
network analysis in this study and previous reports. 
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concentration; and Hi) differentially expressed hypothet- 
ical genes with unknown functions in some of the highly 
correlated modules, especially those in PMT, may play 
important roles in the differential fitness [62]. Based on 
the combined analysis of both modules and hub genes, 
our study indicated that different Prochlorococcus species 
employed some specific metabolic pathways, together 
with general stress responses, to survive in different 
environments. Based our WGCNA analysis results and a 
previous review [63], the representation of responsive 
networks regulated by metals in cyanobacteria was 
schemed in Figure 1. 

Conclusions 

In the study, we have provided evidences that cross- 
species WGCNA analysis could be a powerful tool to in- 
vestigate the molecular mechanisms underlying the 
response of microbes to their environments. Based on 
the shared genes identified between all 4 cyanobacteria 
species evaluated in this study, we have found some 
similar strategies adopted by the cyanoabcteria to deal 
with different stresses. Further investigation of these 
shared genes showing differential expression in various 
species may shed light on their adaption in unique 



environment. In addition, our cross-species analysis has 
found some obvious differences between two metal 
treatments. The extraction of hub genes highly 
connected with other genes in cyanobacteria genomes 
may be valuable in determining essential genes for 
cyanobacteria [64]. Finally, we also noticed many hypo- 
thetical genes with function unknown were grouped 
with important genes in some stress-correlated modules. 
Although they were not discussed in details in this paper 
due to the fact of no much functional information is 
available for them, we believe that further deciphering 
their functions will provide new insights to the stress 
responses in cyanobacteria. 

Methods 

Dataset acquisition 

Microarray datasets of four cyanobacteria species, Pro- 
chlorococcus MED4 (PMM), Prochlorococcus MIT9313 
(PMT), Synechococcus 9311 (SYG) and Synechococcus 
WH 8102 (SYW) [25,31], were downloaded from the 
Gene Expression Omnibus (GEO) [65] (Table 1). The 
experiments were designed for with iron (1 uM) as con- 
trol and without iron as treatment for Prochlorococcus at 
different time (0, 12, 24 and 48 h for PMM and 0, 26, 28 
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Figure 2 Scheme of the WGCNA work flow used in this study. Details of the WGCNA work flow referred in the text and reference (13). 
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and 53 h for PMT based on their growth stages) and 2 h 
of copper treatments with concentrations of 1 \iM 
(SYW) and 0.1 \iM (SYG) [25,31]. The experimental 
details can be found in the original publication [25,31]. 
To compile an extensive set of comparable data, we 
collected all relevant datasets from the two experiments 
and then removed the outliers or samples forming an 
out-group to ensure the results reflecting the real bio- 
logical response to different treatments. Specifically, we 
first performed a hierarchical clustering analysis (ward 
method) using (1-r) as distance for all samples, where r 
is the pair-wise (between any two pairs of samples) 
correlation coefficient for the expression of all genes 
(probesets). Then we compared the tree structure with 
the status of treatment condition. For samples contra- 
dicting with the status of treatment, we excluded them 
from the downstream network analysis. For example, we 
found that in both Prochlorococcus datasets, all samples 
at 0 hour clustered together no matter they were de- 
pleted with iron or had iron, suggesting that their ex- 
pression was irrelevant to iron status. Including them 
will compromise the power of the analysis. Therefore, 
we excluded all samples at 0 hour from our analysis. 
The updated genome annotation for all four species was 
downloaded from NCBI and the Comprehensive Mi- 
crobial Resource (CMR) of TIGR (http://www.tigr.org/ 
CMR) [32-34]. KEGG Pathway and COGs (Cluster of 
Orthologous Groups of proteins) information was obtai- 
ned from the KEGG Pathway Databases (http://www. 
genome.jp/kegg/pathway.html) and the COG Database of 
NCBI (http://www.ncbi.nlm.nih.gov/books/NBK21090/). 

Network construction 

For each cyanobacteria species, we created a trans- 
criptional network from the microarray data, first by 
calculating weighted Pearson correlation matrices cor- 
responding to gene expression, and then by following 
the standard procedure of WGCNA to create the net- 
works [13]. Figure 2 schemes the methodology used in 
this study. Briefly, weighted correlation matrices were 
transformed into matrices of connection strengths using 
a power function [13]. These connection strengths were 
then used to calculate topological overlap (TO), a robust 
and biologically meaningful measurement that encap- 
sulates the similarity of two genes' co-expression re- 
lationships with all other genes in the network [13]. 
Hierarchical clustering based on TO was used to group 
genes with highly similar co-expression relationships 
into modules. From the expression data we followed the 
protocols of WGCNA [60,66] to create within-species 
consensus networks. Gene dendrograms were obtained by 
average linkage hierarchical clustering (Additional file 14: 
Figure S5, A-D), while the color row underneath the den- 
dgram showed the module assignment determined by the 



Dynamic Tree Cut (WGCNA). The network for each mod- 
ule was generated with the minimum spanning tree with 
dissimilarity matrix from WGCNA. To focus on biological 
interpretation, in this study, only genes with functional an- 
notation were included. The modules with r > 0.6, and 
^-value less than 0.05 were extracted for further investiga- 
tion. KEGG pathway analysis was applied to detect the 
common pathways involved among different modules. 
Hub genes were screened by the links (> 8) in the modules 
strongly associated with phenotype (iron depletion or cop- 
per toxicity, based on correlation coefficient r > 0.6). 

To establish cross-species signature genes in this 
analysis, and to define conservation in terms of stress 
response across cyanobacteria species, phenotype-cor- 
related modules and the associated genes in each species 
were extracted. To assess co-expression perseveration be- 
tween the species on a module -by-module basis, we first 
calculated a value of the module membership (MM) - a 
measure of how well each gene correlates with the first 
principal component of gene expression within a module, 
termed as the module eigengene. We then imposed a 
threshold based on MM values (r, p) to make the final 
module assignments. 

Additional files 



Additional file 1: Figure SI. Clustering of the transcriptomic datasets 
under iron and copper treatments in different species. The grouped 
datasets with solid red color in "HighCopper" and "Iron" suggested high 
confident grouping of treated groups from controls. 

Additional file 2: Figure S2. Association between phenotypes and 
identified transcriptional modules in different species. Each of the 
identified transcriptional modules was indicated by different colors, and 
their association with the phenotypes was indicated by p-values. 

Additional file 3: Figure S3. Overlap of the detected responsive 
modules within and between cyanobacteria species. The lines connected 
different cyanobacteria species means there are some genes shared 
between specific modules. 

Additional file 4: Table SI. genes and pathways in SYW-yellow 
module. 

Additional file 5: Table S2. genes and pathways in SYW-Green yellow 
module. 

Additional file 6: Table S3, genes and pathways in SYG-dark green 
module. 

Additional file 7: Table S4. genes and pathways in SYG-brown 
module. 

Additional file 8: Figure S4. Shared responsive genes in cyanobacteria 
species to iron depletion and copper toxicity treatments, respectively, 
and shared responsive genes between treatments in all cyanobacteria 
species. Nine shared genes among all 4 cyanobacteria species were listed 
in Table 3. 

Additional file 9: Table S5. shared genes in iron treatments. The color 
of each gene indicates differential change: red is up-regulated and green 
is down-regulated. 

Additional file 10: Table S6. shared signature genes in copper 
treatments. The color of each gene indicates the differential change, red 
is up-regulated and green is down-regulated. 

Additional file 11: Table S7. Pathways shared in all 4 species. 
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Additional file 12: Table S8. Hub genes with known functions in 
different cyanobacteria species. The same color from KOJD and the 
gene names indicate the hub genes in each species. 

Additional file 13: Table S9. Hub genes with unknown functions in 
different cyanobacteria species. 

Additional file 14: Figure S5. Hierarchical clustering tree using the 
topological overlap dissimilarity. Tree branches have been colored by 
module membership. A, B) PMT, PMM to iron, respectively; C, D) SYW, 
SYG to copper, respectively. Please refer the text for details of the 
analysis. 
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